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Abstract. 

We theoretically examine the rotation of an atomic Bose-Einstein condensate 
in an elliptical trap, both in the absence and presence of a quantized vortex. 
Two methods of introducing the rotating potential arc considered - adiabatically 
increasing the rotation frequency at fixed ellipticity, and adiabatically increasing 
the trap ellipticity at fixed rotation frequency. Extensive simulations of the Gross- 
Pitaevskii equation are employed to map out the points where the condensate 
becomes unstable and ultimately forms a vortex lattice. We highlight the key 
features of having a quantized vortex in the initial condensate. In particular, we 
find that the presence of the vortex causes the instabilities to shift to lower or 
higher rotation frequencies, depending on the direction of the vortex relative to 
the trap rotation. 



PACS numbers: 03.75.Lm,03.75.Kk 
1. Introduction 

Due to their superfluid nature, dilute atomic Bose-Einstein condensates (BECs) are 
constrained to rotate only through the presence of vortices with quantized angular 
momentum [T]. In analogy to the rotating bucket experiment in liquid Helium [2J, 
rotation of a dilute BEC can lead to the nucleation of vortices. The onset of vortex 
nuclcation in a rotating condensate is non-trivial and has been the subject of much 
experimental H ED E3 EH and theoretical Dfl El USB ttH 021 031 031 0S1 OSl HZI 
[TBI [THl HOI [H] investigation. One method of spinning the condensate is to use a 
rotating elliptical harmonic trap, with the presence of vortices predicted to become 
energetically feasible for rotation frequencies exceeding f2 ~ 0.3wj_ [5], where uj± is 
the average trap frequency in the rotating plane. Experimentally, however, vortices 
are formed only when the trap is rotated at much higher frequencies O 01 [5] . This 
anomaly exists because it is only at these higher frequencies that a dynamical route 
to vortex nucleation appears. 

In the rotating frame, the static condensate solutions in an elliptical trap are 
a family of quadrupole solutions with elliptical density profiles. These solutions, 
and their regimes of stability, can be conveniently approximated using the classical 
hydrodynamic equation in the rotating frame [12[ 113] . Experimentally, these solutions 
can be accessed by two methods: Procedure I involves increasing the rotation 
frequency from zero while maintaining a fixed elliptical trap, while Procedure II 
involves increasing the trap ellipticity from zero while rotating at a fixed frequency. 
Providing these procedures are performed adiabatically, the condensate will access the 
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stationary states in the rotating frame. Indeed, this has been observed experimentally 
[7j- However, depending on which procedure is employed, a critical rotation 
frcqucncy/cllipticity is reached beyond which these vortex-free solutions become 
unstable and lead to the formation of a vortex lattice. The onset of the instability has 
been mapped out theoretically based on classical hydrodynamics [12l HSJ [20] and is 
in excellent agreement with experimental results [5] [7]. Furthermore, time-dependent 
simulations of the highly-successful Gross-Pitaevskii equation [22] have elucidated the 
full dynamics of this process [H] [HI [16] [17l HH [19] [20], including the non-trivial 
transition from instability to vortex lattice. If the rotating elliptical trap is introduced 
non-adiabatically, instead, the condensate undergoes shape oscillations which can be 
unstable and lead to vortex lattice formation. These non-adiabatic dynamics have 
been observed experimentally [7] and elucidated theoretically (14 ] [18 ] [19] . 

A natural extension of these studies is to consider how the presence of a quantized 
vortex in the initial condensate affects the dynamics under rotation. To date this has 
been considered by a handful of theoretical studies. The effect of a vortex on the 
collective modes of a condensate has been analysed (23] [24] [25] [26] and predicted to 
induce an upwards shift in the mode frequencies. Tsubota et al. [14] showed that a 
rotating condensate containing a vortex undergoes similar dynamics to the vortex-free 
case, although the final configuration of vortices depends on the initial state. Such 
hysteresis effects, depending on the initial angular momentum state of the condensate, 
have been studied in more detail by Jackson and Barenghi [21] . Crucially, these 
theoretical studies consider only the case of a vortex whose flow is concurrent with 
the trap rotation. In this paper we consider the rotation of a condensate containing a 
vortex whose flow is either concurrent or against the trap rotation. We examine the 
adiabatic introduction of the rotating elliptical trap, cither by increasing the rotation 
frequency for a fixed trap ellipticity (Procedure I) or by increasing the trap cllipticity 
at fixed trap rotation frequency (Procedure II). We map out the onset of instability 
that triggers vortex lattice formation for both of these procedures, and compare to 
the non-vortex case. 

In section[5]we outline our theoretical approaches, namely, the full time-dependent 
simulation of the condensate using the Gross-Pitaevskii equation and an analysis of 
the classical hydrodynamic solutions in the rotating frame. In section [3] we describe 
the instabilities of a vortex-free condensate that arise from introducing a rotating 
elliptical trap, and present new results. In particular, we make a thorough comparison 
between the analytic hydrodynamic predictions and the computational simulations of 
the Gross-Pitaevskii equation. In section 3] we examine how the presence of a vortex 
changes these dynamics, and, in particular, contrast the cases where the vortex flow 
is concurrent and against the rotating trap. Finally, in section [5] we present the 
conclusions of our investigation. 

2. Theory 

Assuming ultra-cold temperature and weak atomic interactions, the vast majority 
of the bosons in the system are in the Bose-condensed state. Thermal and 
quantum effects can become negligible and the condensate can be parameterised by 
a macroscopic "wavefunction" ip(r,t). Moreover, this wavefunction is governed by 
the mean- field Gross-Pitaevskii equation (GPE). In the frame rotating at frequency 
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Q about the z-axis this equation has the form, 



ift 



2m 



V 2 + y cxt (r,t)+.g|^(r, £)| 2 - Oi, 



V'M), (1) 



where T4 x t(r,£) is the external confinement and m is the atomic mass. The nonlinear 
coefficient g = Anh 2 Na/m characterises the atomic interactions, where N is the 
number of atoms and a is the s-wave scattering length. The term containing the 
angular momentum operator L z accounts for frame rotation about the z-axis. The 
GPE represents the back-bone of theoretical studies into BECs and has been shown 
to give an excellent description of many static and dynamical effects, including vortices 

mm. 



2.1. Hydrodynamic approach 

It is useful to consider the hydrodynamic interpretation of the wavefunction ijj(r, t) = 
\J p(r, t) exp[i</>(r, t)], where p(r,t) is the atomic density and </>(r, i) is a phase factor 
which defines a fluid velocity via v(r, t) = (% J 'm)V 4>{r , t) . Substitution of these 
relations into the GPE and equating real and imaginary parts leads to the classical 
hydrodynamic equations, 

^+V-[p(v-fixr)]=0, (2) 

m ~dl + V (^"^ ' V + ^ + pg ~ mv ' ^ X = °' 

Equation @ is a continuity equation while equation (J3|) has the form of a Euler 
equation for an irrotational fluid. In deriving equation we have employed the 
Thomas-Fermi approximation by neglecting the non-trivial 'quantum pressure' term 

fi^ nana ess. 

The confining potential V cx t (r 5 1) is typically formed by magnetic fields and has 
the harmonic form, 

V ext (r,t) = \m{u? x x 2 + u 2 y y 2 + uj 2 z z 2 ), (4) 

where u> x , oj v and uj z are the harmonic trap frequencies in the three cartesian 
dimensions. In the experiments of [3l [5l [Z] the trap is approximately elliptical in the 
rotating (x—y) plane. Therefore we employ an cllipticity parameter e to define the x, y 
trap frequencies to be u x = \/l — e u>± and u> y = y/l + e lo±, where <J\_ = {lo 2 +u> 2 )/2 
is the geometric mean of the frequencies. 

Rotation of the elliptical trap tends to excite a quadrupolar mode in the 
condensate. Following Recati et al. [T^, we assume an irrotational quadrupolar flow 
in the condensate defined by the velocity field [53] > 

v = aV(xy), (5) 

where a is the velocity field amplitude. It is important to note that a is also 
proportional to the cllipticity of the condensate density profile. Following the analytic 
procedure in references [T^l EH [20] , the classical hydrodynamics equations leads to a 
density profile given by, 

1( 1 /~2 2,~2 2, 2 2\\ / a \ 

p= -\u- -m{w x x +w y y +w s z)\, (6) 
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Figure 1. The size of the largest real part of eigenvalues found using equation 
J5). The size of the eigenvalue is indicative of the speed at which the condensate 
breaks down upon entering an unstable region. The inset shows the same plot but 
at much lower amplitude to highlight the unstable regions with small eigenvalues 
which are not visible in the main plot. 



for p > 0, and p = elsewhere. Furthermore, this leads to the equation, 

a 3 + (1 - 2fl 2 )a - efl = 0, (7) 

which defines the stationary condensate solutions (defined by the quadrupolar flow 
amplitude a) in a harmonic trap with ellipticity e rotating at frequency Jl. Throughout 
this paper we will consider the relevant parameter space (e — Q) which defines these 
condensate solutions. 

Although equation ([7]) defines static solutions, they are not necessarily stable 
solutions. To examine their dynamical stability, we consider small perturbations 
in density Sp and phase 5<fi to the static solutions [131, ED] ■ Taking the variational 
derivatives of equations ([2]) and (|3|) leads to the time-evolution equations 

g/fi 

v c - V 



d_ 

dt 



' Scf> ' 
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(8) 



v c • V 

where v c = v — $7 x r is the velocity field in the rotating frame. A polynomial ansatz is 
taken for the perturbations Sp and S(p, and inserted into equation ([5]). As a technical 
note, our polynomial ansatz consists of terms of the form x l y^ where i,j > and 
i+j < 7. If the resulting eigenvalues are purely imaginary (or contain a negative real 
part) then the perturbations are stable oscillations (or decay over time). However, 
if any of the eigenvalues have a real positive component, the perturbations will grow 
exponentially over time. Such solutions are dynamically unstable and lead to vortex 
lattice formation [T31 HO] . 

Equation @ defines regions in (e — £1) space that arc unstable as shown in figure 
[1] When the ellipticity and rotation frequency of the trap enters this region the 
instability leads to disruption of the condensate, vortex nuclcation and ultimately 
lattice formation. The speed at which the instability evolves will increase with the 
size of the real part of the eigenvalues. The size of the largest real component of the 
eigenvalues are mapped out in figure [T] This shows that some regions of instability will 
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take significantly longer to produce vortices than others. This effect will also studied 
in detail in section [31 

Another source of instability to consider is the center of mass instability which 
occurs when the rotation frequency becomes close to the trapping frequency [12] . This 
instability is equivalent to the classical case of a point particle in a rotating elliptical 
potential. For rotation frequencies f2 lying between uj x and uj y the oscillations of 
the trapping potential couple to the oscillation frequencies of the particle and it is 
ejected from the trap. This has been experimentally observed in the explosion of the 
condensate [57]. However, for sufficiently strong repulsive interactions, explosion is 
prevented and the condensate remains intact, although its centre of mass deviates 
from the trap centre [T9J, [27] . 

2. 2. Numerical solution of the GPE 

Although the hydrodynamic approach outlined above gives valuable insight into the 
rotating condensate solutions, it has limitations. For example, it assumes the Thomas- 
Fermi limit, it can only predict the stability of modes with known flow patterns, and 
it cannot give information about how the instability manifests itself or the final state 
of the system. A more thorough method, albeit more intensive, is to explicitly solve 
the full Gross-Pitaevskii equation. Importantly, this method allows us to simulate the 
time dynamics of the condensate from the initial state all the way to its final state. 

We employ the Crank-Nicholson method to numerically propagate the Gross- 
Pitavcskii equation in time. The dynamics predicted by equations ([5]) and (J8]) are 
independent of z 120] so we are justified in conducting the simulations in two- 
dimensions. Condensate solutions are found by the imaginary-time technique. Under 
the substitution (t — ► —it) and using an appropriate initial guess, the GPE converges 
to the lowest energy solution of the system, provided it is stable [21] ■ 

In presenting our numerical results, we employ the healing length £ = h/ y/rrmog 
as the unit of length and (£/c) as the unit of time, where no is the peak condensate 
density of the initial (non-rotating) solution and c = y nog/m is the Bogoliubov speed 
of sound. The chemical potential fx = n$g at t = is the unit of energy. We consider 
a condensate with /i = 7h~uj and a Thomas-Fermi radius of i?TF = 10£. Assuming 
87 Rb atoms, our units of space and time typically correspond to £ ~ 0.2 fim and 
(C/c) ~ 10- 4 s. 

Once the solutions are found, they are propagated with the GPE in real time. We 
follow the experimental procedures as closely as possible [U [5] [7] . We first obtain the 
initial solution with either f2 = or e = and then introduce the rotating elliptical trap 
adiabatically by following Procedure I or Procedure II. The condition of adiabaticity 
is important because it means that the evolving condensate will access the rotating 
frame solutions. The non-adiabatic introduction of the rotating elliptical trap leads, 
instead, to oscillations between a non-elliptical and elliptical state [15 ] \19 \ . 

Recent numerical studies [TH] [TH] [5D] have highlighted the necessity of rotational 
symmetry breaking for successful and realistic vortex nuclcation. For example, if the 
two-fold rotational symmetry of the system is preserved then vortices must enter in 
opposing pairs, which is highly unfavourable in energetic terms. Consequently, at 
the start of each simulation of the GPE we shift the trap centre by a distance of 
0.1£. This breaks the two-fold rotational symmetry of the system and enables realistic 
vortex nuclcation. 
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Figure 2. (a) Hydrodynamic quadrupolc condensate solutions a of the rotating 
condensate as a function of rotation frequency according to equation J7J for e = 
(solid line), e = 0.02 (dotted line) and e = 0.1 (dashed line), (b) The position 
of the bifurcation point £lb( € ) (solid line) and the unstable region (grey area) of 
the hydrodynamic condensate solutions are indicated. The crossing frequency f2x 
is shown. In both (a) and (b), the routes for introducing the rotating elliptical 
trap, Procedure I and II, are illustrated. For Procedure II, rotation in the regime 
f2 < fl^ accesses the upper branch of the solutions, while rotation in the regime 
$1 > accesses the lower branch. 



2.3. Energy of a rotating condensate 



In the laboratory frame the energy is not conserved due to the presence of the time- 
dependent potential. The energy in the rotating frame is conserved (providing no 
dissipation is present) and is given by, 
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The terms in the integral represent the kinetic energy Ek, potential energy Ep and 
interaction energy E\. The final term represents the rotational energy, where the 
average angular momentum per particle about the z-axis, (L z ), is given by, 
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3. Rotating a vortex-free condensate 

The instability of a vortex-free condensate under the adiabatic introduction of 
a rotating elliptical trap has been mapped out by experiment [5J [7] and theory 
[T2l [TBI, UM |2"U] . We will outline the key features of these studies, before presenting 
additional results. We generally consider the frequency regime < u)±. Above 
this frequency, centrifugal forces overcome the trap potential and the condensate can 
undergo explosion or centre-of-mass oscillations [12j [19j [27] . 

Figure (2]^a) shows the static solutions in (a — f2) space according to equation j7]). 
It should be recalled that a defines the condensate solution via equation ([5j and is 
proportional to the ellipticity of the BEC density profile. For an axially-symmetric 
trap e = (solid line), the a = solution exists throughout while for O > uj±/^/2 two 
additional solutions bifurcate symmetrically to finite negative and positive values of 
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a. We denote the frequency of this bifurcation point to be 0°. For finite e [dashed 
and dotted lines in figure 0Ja)] there exist two distinct branches of solutions. The 
upper branch of a is positive, single-valued and occurs for all rotation frequencies. 
The lower branch is negative, double- valued, and only exists for > Ob(e)- The 
bifurcation frequency fib(e) is now a strict function of ellipticity and increases with e 
as plotted in (e — ft) in figure Elb) (solid line). 

In addition to the static solutions, we also need to consider their dynamical 
instability according to equation (JSj) . This predicts that the condensate solutions are 
dynamically unstable within a region of (e — £1) space, as shown by the shaded area in 
figure Hlb). As e is increased, the critical rotation frequency for dynamical instability 
becomes lower than fib(e)- We define the frequency at which these instabilities cross 
to be the crossing frequency Q x = 0.765a; j_. 

We are now in a position to predict the behaviour of the condensate as a rotating 
elliptical trap is introduced, with the results presented below. 

3.1. Procedure I 

From the hydrodynamic point of view, the trap initially has fixed ellipticity, e.g. 
e = 0.1 (dashed line in figure^a)). As the rotation frequency fl is increased from zero, 
the condensate follows the upper branch solution, becoming increasingly elongated as 
it moves to higher values of a. As the BEC moves horizontally across figure [DJb), it 
becomes dynamically unstable when it reaches the unstable region (shaded region of 
figure Hb)). 

We have performed numerical simulations of the condensate while Q is 
adiabatically introduced. An example is presented in figurc[3]for a fixed trap ellipticity 
of e = 0.05. The rotation frequency f2(i) is ramped up linearly from to 0.8lo± over 
a time of 2000(£/c) and then maintained at a fixed value, as shown in figure [3Ja). 

Figure [3je) shows snapshots of the condensate density during its evolution, while 
figure [3]Jf) shows the corresponding plots of the condensate phase. The condensate is 
initially circular in shape [figure [3je) (i)] with a constant phase [figure [3ff)(i)]. As the 
rotation frequency is increased from zero, the condensate becomes elongated [figure 
[3fc(ii)] as it traces out upper branch of the hydrodynamic solutions. The irrotational 
velocity field of equation ([5]) correspond to a phase profile <fi(x, y) oc xy which is clearly 
evident in figure [3f(ii). On the outskirts of this irrotational phase pattern we see 2tt 
phase singularities. These are termed 'ghost' vortices because they cannot be seen in 
the density profile [14] . However, at some point the 'ghost' vortices penetrate the bulk 
of the condensate such that it begins to deviate from the smooth quadrupole solution 
[figure [He-f)(iii)], marking the onset of the dynamical instability. This instability 
disrupts the condensate into a turbulent state [figure [3je) (iv)] , which subsequently 
relaxes into a vortex lattice [figure Efe)(v)] via vortex-sound interactions [18] . 

At early times, the z-component angular momentum per particle (L z ) is zero, 
as shown in figure [3Jb). However, when the instability kicks in angular momentum 
becomes rapidly driven into the system by the rotating trap, up to some maximum 
value. Figure [3jc) shows the evolution of the kinetic energy (solid line), interaction 
energy (dashed line) and potential energy (dotted line) during the dynamics. We have 
normalised these positive components of the energy by the laboratory frame energy 
El = (Ek + Ep + Ei). 

The condensate instability is associated with a sharp increase in the kinetic energy 
as the trap drives the condensate into an energetic turbulent state, and a decrease in 
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Figure 3. (a) The ramping up of the trap rotation frequency over time, (b) 
The average z-componcnt angular momentum per particle during the simulation, 
(c) The evolution of the condensate kinetic energy (solid line), interaction energy 
(dashed line) and potential energy (dotted line) under the rotation frequency 
ramp shown in (a) and for a fixed trap ellipticity e = 0.05. (d) Snapshots of the 
condensate at (i) t = 0, (ii) t = 1775, (iii) t = 2320, (iv) t = 3670, and (v) 
t=15000(£/c). Each image is of dimension 50£ by 50£, although the simulations 
were conducted in a numerical box of size 100£ by 100£. (e) Corresponding 
snapshots of the condensate phase. 



the interaction energy as the condensate becomes more diffuse. 

Following the instability, the dynamics are explicitly time dependent and progress 
regardless of whether the trap frequency continues to be ramped up or is stopped. As 
such, we employ this characteristic time dependent deviation of the energies to define 
the point of condensate instability. In figure for example, the instability kicks in at 
t « 1800(£/c), corresponding to a rotation frequency of O « 0.73uj±. With a very slow 
ramping rate this point can be pinpointed to very high accuracy and is more exact 
than in previous studies where the point of instability was determined 'by eye' [20] . 

Using the GPE we have mapped out the rotation frequency at which the 
condensate becomes unstable for a range of trap ellipticitics and for two different 
ramping rates, with the results plotted in figure [DJa) (circles and squares). We find 
that the ramping rate can significantly affect the point at which the condensate breaks 
down. For e < 0.05 we employ a ramping rate of dil/dt = 2 x 10 _4 c^ (squares). 
Here the GPE shows that the condensate becomes unstable in the frequency range 
0.7<x>j_ ^ Q ^ 0.8u>±, in good agreement with the boundary of the main dynamically 
unstable region (shaded region). For e > 0.04 we employ a slower ramping rate of 
dfl/dt = 4x lO -6 ^ (circles). Here the GPE simulations show that the condensate 
instability sets in at much lower rotation frequencies (0.55u)±^^1^0.6uj±) that are 
consistent with the narrow band of dynamical instability. As discussed in section 
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P| this narrow band of dynamical instability features small eigenvalues, and so the 
dynamical instability is slow to evolve. The simulations show that, providing the 
ramping rate is high enough, the condensate can pass through the narrow region of 
dynamical instability with negligible effect. Here we have not presented the slow 
ramping rate for cllipticitics e < 0.04 since the point of instability becomes less clear. 
Due to the weakened effect of the narrow instability region on the condensate (reduced 
eigenvalues and narrower width), the condensate passes through the narrow region of 
instability before becoming fully unstable. However it undergoes a partial instability 
that distorts the density slightly from its initially smooth profile. Upon emerging from 
this unstable region the condensate slowly relaxes to its steady state solution. When 
it reaches the second instability region the condensate breaks down fully. For slower 
ramping rates the condensate will break down fully within the narrow unstable region 
for e < 0.04. For the same reason our data for the high ramping rate is limited to 
cllipticities e < 0.05. 

Madison et al. [7] adiabatically increased the rotation frequency at a fixed 
cllipticity of e w 0.028 and observed the onset of condensate disruption at rotation 
frequency SI ss 0.75wj_. This point, indicated on figure [4]J a) by the cross, is in good 
agreement with the theoretical predictions, particularly the GPE results. 

3.2. Procedure II 

In this method the trap cllipticity is increased adiabatically from zero while the 
rotation frequency is kept fixed. Procedure II has been employed experimentally 
to generate vortices by Hodby et al. [5] and theoretically analysed using the 
hydrodynamic approach [12j [13] and GPE simulations [16l [20] . In these latter 
studies HO], the GPE results were found to be in very good agreement with 
the hydrodynamic predictions and experimental results [5]- In particular, reference 
[20] elucidated the condensate instabilities that arise from Procedure II. Up to some 
critical ellipticity, the condensate accesses stable quadrupole solutions. However, once 
the critical ellipticity is reached the condensate becomes disrupted. There are three 
distinct regimes of instability depending on SI: 

• Ripple instability fi < < uj±/\/2: The condensate follows the upper branch 
solutions to increasing values of a, as indicated in figure [2] by the arrow labelled 
II(tt < Jig). However, at some critical ellipticity e c the solutions become 
dynamically unstable according to equation ©. GPE simulations reveal that 
the dynamical instability manifests itself in the appearance and growth of density 
ripples, which ultimately lead to the complete disruption of the condensate. 

• Interbranch instability Slg < SI < Sl x : The condensate follows the lower branch 
solutions to negative values of a with increasing magnitude, as indicated by the 
arrow labelled //(SI > Sl°)- Eventually the lower branch ceases to be a solution. 
Since the upper branch solutions are stable the condensate tries to deform to 
the upper branch solutions, generating large unstable shape oscillations, which 
disrupt the condensate. 

• Catastrophic instability SI > Sl x : As with the interbranch instability the 
condensate follows the lower branch solutions (route //(SI > Sl£) in figure^. 
However, when the lower branch solution disappears, the upper branch is also 
unstable, and so no stable solutions exist. A large and sudden contortion of the 
condensate occurs, leading to a highly disrupted state. 
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Figure 4. Phase diagrams in (e — Q) space showing the instability of a rotating 
elliptical condensate for (a) Procedure I (fixed e and increasing Q). Squares 
represent the onset of instability according to the GPE simulations for a fast 
ramping rate of dCl/dt = 2 X 10 -4 c^, while circles represent a slow ramping 
rate of dQ/dt = 4 X I0~ e uj 2 ± . The cross is the experimental data from [7]. (b) 
Procedure II (fixed Q and increasing e), the frequency of the bifurcation instability 
is givenSlb( e ) (solid line). The crosses represent the experimental data of [5]. The 
dynamical instability (grey area) is shown in both plots. 



Each of these regimes leads to a turbulent state of vortices and sound waves, which 
relaxes into a vortex lattice via vortex-sound interactions [T5] . 

Figure[4th) maps out the point of the condensate instability in (e — f2) space under 
Procedure II. This figure is similar to figure 1(b) of reference [20], but employs an 
improved definition of the point of instability. Procedure II involves moving vertically 
upwards in (e — O) space from the e = axis. Based on the hydrodynamic solutions 
and their dynamical stability, we expect that for Q < fix the condensate will become 
unstable when it reaches the bifurcation point (solid line) while for fl > Six the 
condensate will become unstable when it reaches the dynamically unstable region 
(shaded region). 

We have performed GPE simulations of the condensate while the ellipticity is 
ramped up linearly, over a range of rotation frequencies. To ensure adiabaticity we 
employ a slow rate of change of ellipticity, de/dt = 10~ 4 cjx- When the state of the 
condensate becomes explicitly time dependent, we have isolated the critical ellipticity. 
These GPE results, shown by the circles in figure HJb), are in very good agreement 
with the hydrodynamic predictions, despite the assumptions of the hydrodynamic 
model. The experiment of Hodby et al. [5] performed Procedure II to investigate the 
instability of the condensate. The experimental results, shown by crosses in figurc|D(b) , 
are well-described by the both theoretical methods, particularly the GPE method. 

4. Rotating a condensate containing a vortex 

We now consider the effect of rotating a condensate that already contains a singly- 
quantized vortex at its center. The hydrodynamic approach outlined in section 
12.11 assumes a node-less Thomas-Fermi density profile and a quadrupolar velocity 
field. These approximations are no longer valid in the presence of a vortex, which 
creates a non-Thomas-Fermi node in the density and modifies the velocity field of the 
condensate. For this reason we proceed by solving the full GPE numerically, and make 
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Figure 5. Instability of a condensate containing a vortex following Procedure 
I with a fixed trap ellipticity of e = 0.085. (a) z-componcnt angular momentum 
per particle (L z ) as the rotation frequency Q is introduced in the presence of a 
vortex concurrent with (dashed line) and against (dotted line) the trap rotation. 
The results for a vortex- free condensate are also shown (grey line), (b) Density 
and (c) phase snapshots of the dynamics for the case of a vortex against the trap 
rotation. These snapshots correspond to rotation frequencies of (i) Q = 0.008wx > 
(ii) Q = 0.46w ± , (iii) U = 0.47oj x and (iv) fi = 0.54w x . 



comparisons to the non- vortex case which is understood analytically. In order to form 
the vortex state, we enforce a vortex phase profile 4>(x,y) = tan _1 (y/x) during the 
imaginary time propagation method. A typical initial state is shown in figure [5jb) (i) 
and (c)(i). 

4-1- Procedure I 

Using Procedure I, we rotate a condensate containing a vortex whose circulation is 
either concurrent with or against the trap rotation. Figure [5] presents the typical 
evolution of the condensate for a fixed trap ellipticity of e = 0.085. For the vortex- 
free condensate (grey line in figure [5]), the growth of the angular momentum, which 
indicates the onset of instability, occurs at « 0.6wj_. When the vortex is concurrent 
with the trap rotation (dashed line), the angular momentum per particle is L z = h 
at early times. The onset of instability occurs at a higher frequency, f2 ss 0.65wj_. 
The instability progresses in a similar manner to the vortex-free case. 'Ghost' vortices 
develop on the condensate edge with the same flow direction as the initial vortex and 
the trap rotation. Eventually the state of the condensate becomes unstable. Once 
the instability starts, the growth of angular momentum occurs at a similar rate to the 
vortex-free case. 

For the case of a vortex flowing against the trap rotation, the initial angular 
momentum is L z = —%. Once the trap rotation frequency is increased from zero, the 
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Figure 6. Instability of a condensate containing a vortex following Procedure I. 
Phase diagram in (e — Q) space showing the onset of condensate instability under 
Procedure I for a vortex concurrent with (squares) and against (circles) the trap 
rotation. The instability points for a vortex-free condensate are also presented 
(grey stars). 



presence of this oppositely-rotating vortex is highly unfavourable in the system. At 
a relatively low frequency, which is f2 « 0.4cjj_ for the example shown in figure [5ja), 
this configuration becomes dynamically unstable. Due to the highly unfavourable 
configuration, the instability progresses rapidly. 'Ghost' vortices appear on the 
condensate edge, whose circulation is concurrent with the trap rotation and in the 
opposite direction to the central vortex [see the density and phase images in figures 
[5jt>) (ii) and (c)(ii)]. Their presence increases the angular momentum. Once the 
instability is reached, one or more ghost vortices become driven towards the condensate 
center [figures [5^a) (iii) and (b)(iii)]. For the example shown in figure [5] the opposing 
vortex becomes ejected from the condensate and the system settles into a configuration 
containing one concurrent vortex [figures [SJa)(iv) and (b)(iv)] and positive angular 
momentum of the order of unity. In other words, this initial phase of the instability 
has reversed the vortex configuration to the more stable, concurrent configuration. 
The system then progresses in this concurrent vortex state, remaining approximately 
stable until it reaches a second, higher critical frequency at f2 ~ 0.7lu±. Angular 
momentum and multiple vortices becomes driven into the condensate and ultimately 
form a vortex lattice. 

Using the GPE we have measured the onset of instability of Procedure I over 
a range of trap cllipticities, with the results presented in figure [H We consider the 
case where the vortex flow is concurrent with (squares) and against (stars) the trap 
rotation. The results from section 13.11 for a vortex-free condensate are shown for 
comparison (grey circles). The behaviour of the critical rotation frequencies for each 
case is non-trivial. Crucially, however, the presence of the vortex shifts the critical 
frequencies. For a vortex concurrent with the rotation, the critical frequencies become 
shifted to higher values, while for a vortex against the rotation, the critical frequencies 
are shifted to lower frequencies. 

In the absence of a vortex, the quadrupolc mode of a condensate is predicted 
to become unstable at a rotation frequency ft = uj±/^/2 [28] . This assumes an axi- 
symmctric system (e — 0) in the Thomas-Fermi limit. In a similar manner, several 
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Figure 7. (a) Phase diagram in (e — f2) space showing the onset of condensate 
instability under Procedure II for a vortex concurrent with (squares) and against 
(stars) the trap rotation. The instability points for a vortex-free condensate are 
also presented (grey circles), (b) The centre-of-mass instability in the region 
f2 ~ uj x ,yi with the analytic prediction of e = ±(1 — O 2 ) shown (solid line) 1121 . 



studies have theoretically considered the quadrupole mode of a condensate containing 
a vortex [531 [Ml HS1 HI] ■ It is predicted that the presence of the vortex shifts the 
critical frequency by Af2 = (77kJj_/8/ti). For the condensate parameters employed in 
this study (fi = 7hu)±), this shift is Af2 = 0.125o;_l, which is in good agreement with 
the observed shifts in figure [5] 

4-2. Procedure II 

We now consider a condensate containing a vortex and following Procedure II. The 
critical cllipticitics at which the condensate becomes unstable are mapped out in figure 
[TJa) for a vortex concurrent with (squares) and against the trap rotation (stars). The 
results from section 13.21 are shown for comparison (grey circles) . Figure \T[b) is the 
region of the plot centered around the center of mass instability (where fl « ui±). The 
solid line shows the theoretical position of the center of mass instability. All the cases 
(regardless of vortex configuration) break down at the same point, as predicted by the 
theory (as this phenomenon is classical, it depends only on the center of mass (x) and 
not on the quantum mechanical distribution of the wavefunction) 

Again, the condensate instability is shifted to higher frequencies by the presence 
of a concurrent vortex and lower frequencies by the presence of an oppositely rotating 
vortex. 

Although the presence of a vortex and its configuration within the BEC has a large 
effect on the stability of the system, once the BEC has become unstable its behaviour 
is similar in all cases. The kinetic energy's contribution to the total energy increases 
substantially whilst the contribution from the interaction energy becomes much less 
important. In all cases the instability allows vortices to enter the condensate, which 
ultimately settle into a lattice. 

5. Conclusions 

In this paper we have analysed and mapped out the instability of a condensate to 
a rotating elliptical trap, both in the absence and presence of a vortex in the initial 
condensate. We follow two methods of introducing the driving potential - Procedure I 
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Figure 8. Angular momentum per particle (L z ) as Procedure II is followed in 
the presence of a vortex concurrent with (dashed line) and against (dotted line) 
the trap rotation. The results for a vortex-free condensate are also shown (grey 
line). The rotation frequency is (a) Q = 0.6uj± an d (b) ^ = 0.85wx- 



involves increasing the rotation frequency at fixed trap cllipticity, while Procedure II 
involves increasing the trap ellipticity at fixed rotation frequency. These procedures 
are performed adiabatically so that the stable condensate accesses the quadrupolar 
static solutions in the rotating frame. In the vortex-free case, these can be analysed 
within the Thomas- Fermi limit using the classic hydrodynamic equations. At a 
critical rotation frequency/ellipticity, the solutions become unstable. We map out this 
point of instability using time-dependent simulations of the Gross-Pitaevskii equation. 
Although the non-vortex case has been examined in detail elsewhere, we present 
some additional data and demonstrate the good agreement with the hydrodynamic 
predictions and experimental data. 

In the presence of a vortex, the points of instability become shifted. For a vortex 
which is concurrent with the trap rotation, the instability is shifted to higher rotation 
frequencies, while for a vortex which is against the trap rotation, the instability is 
shifted to lower rotation frequencies. The shift is of the order of 0.1lu±, which is in 
good agreement with the shift of the quadrupolc mode that is predicted to occur in 
the presence of a vortex. 
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